####### post #######

setGeneric("post", 
           function(model,x1name=NULL,x1vals=NULL,x2name=NULL,x2vals=NULL,holds=NULL,seed=NULL,
                    n.sims=1000,cut=NULL,quantiles=c(.025,.975),did=NULL,weights=NULL, digits=2){
             standardGeneric("post")
           }
)

setClassUnion("arrayORNULL", c("array","NULL"))
setClassUnion("listORcharacter", c("list","character"))

setClass("post", 
         slots = c(est = "array", 
                   did = "arrayORNULL", 
                   sims = "array", 
                   model = "character", 
                   link = "listORcharacter", 
                   quantiles = "numeric", 
                   call = "call")
)


post.glm <- function(model,x1name=NULL,x1vals=NULL,x2name=NULL,x2vals=NULL,holds=NULL,seed=NULL,
                     n.sims=1000,cut=NULL,quantiles=c(.025,.975),did=NULL,weights=NULL, digits=2){
  
  call <- match.call()
  
  set.seed(seed)
  sims <- postSim(model, n.sims=n.sims)
  
  if (family(model)[2]=="identity"){link <- identity} 
  else if (family(model)[2]=="probit"){link <- pnorm}
  else if (family(model)[2]=="logit"){link <- plogis}
  else if (family(model)[2]=="log"){link <- exp}
  else if (family(model)[2]=="cloglog"){link <- function(x){1-exp(-exp(x))}}
  else {stop("Link function is not supported")}
  
  if (is.null(weights)){wi <- c(rep(1, length(model$model[,1])))} else{wi <- weights}
  n.obs <- length(model.matrix(model)[,1])
  k <- length(model.matrix(model)[1,])
  n.q <- length(quantiles)
  
  if (is.null(x1name)){
    X <- array(NA, c(n.obs,k))
    newdata <- data.frame(model$model)
    if (!is.null(holds)){
      for (i in 1:length(holds)){
        newdata[ ,names(holds)[i]] <- as.numeric(holds[i])
      }
    }
    X <- t(model.matrix(lm(formula(model), data=newdata)))
    l1 <- array(NA, c(nrow(sims@coef),1))
    l1[,1] <- apply(link(sims@coef %*% X), 1, function(x) weighted.mean(x, wi))
    l2 <- array(NA, c(1,n.q+1))
    l2[1,1] <- weighted.mean(link(coef(model) %*% X), wi)
    l2[1,2:(n.q+1)] <- quantile(l1, probs=quantiles)
    colnames(l2) <- c("est",quantiles)
    
    ans <- new("post", 
               est=round(l2, digits=digits),
               did=NULL,
               sims=l1, 
               model=class(model), 
               link=family(model)[2], 
               quantiles=quantiles, 
               call=call)
    return(ans)
  }
  
  else if (is.null(x2name)){
    
    n.x1 <- length(x1vals)
    X <- array(NA, c(n.obs,k,n.x1))
    
    for (i in 1:(n.x1)){
      newdata <- data.frame(model$model)
      if (!is.null(holds)){
        for (j in 1:length(holds)){
          newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
        }
      }
      newdata[ ,x1name] <- x1vals[i]
      X[ , ,i] <- model.matrix(lm(formula(model), data=newdata))
    }
    
    X <- aperm(X, c(2,1,3))
    l1 <- apply(apply(X, c(2,3), function(x) link(sims@coef %*% x)), c(1,3), function(x) weighted.mean(x, wi))
    l2 <- array(NA, c(n.x1+1,n.q+1))
    l2[1:n.x1,1] <- apply(apply(X, c(2,3), function(x) link(coef(model) %*% x)), 2, function(x) weighted.mean(x, wi))
    l2[1:n.x1,2:(n.q+1)] <- aperm(apply(l1, 2, function(x) quantile(x, probs=quantiles)))
    l2[nrow(l2),1] <- l2[n.x1,1] - l2[1,1]
    l2[nrow(l2),2:(n.q+1)] <- quantile(l1[ ,n.x1] - l1[ ,1], probs=quantiles)
    rownames(l2) <- c(paste(c(rep(paste(x1name,"="),n.x1),
                              paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),
                            c(x1vals,"")))  
    colnames(l2) <- c("est",quantiles)
    
    ans <- new("post", 
               est=round(l2, digits=digits), 
               did=NULL, 
               sims=l1, 
               model=class(model), 
               link=family(model)[2], 
               quantiles=quantiles, 
               call=call)
    return(ans)
  } 
  
  else{
    
    n.x1 <- length(x1vals)
    n.x2 <- length(x2vals)
    X <- array(NA, c(n.obs,k,n.x1,n.x2))
    
    for (j in 1:n.x2){
      for (i in 1:n.x1){
        newdata <- data.frame(model$model)
        if (!is.null(holds)){
          for (k in 1:length(holds)){
            newdata[ ,names(holds)[k]] <- as.numeric(holds[k])
          }
        }
        newdata[ ,x1name] <- x1vals[i]
        newdata[ ,x2name] <- x2vals[j]
        X[ , ,i,j] <- model.matrix(lm(formula(model), data=newdata))
      }
    }
    
    X <- aperm(X, c(2,1,3,4))
    l1 <- apply(apply(X, c(2,3,4), function(x) link(sims@coef %*% x)), c(1,3,4), function(x) weighted.mean(x, wi))
    l2 <- array(NA, c(n.x1+1,n.q+1,n.x2))
    l2[1:n.x1,1,1:n.x2] <- apply(apply(X, c(2,3,4), function(x) link(coef(model) %*% x)), c(2,3), function(x) weighted.mean(x, wi))
    l2[1:n.x1,2:(n.q+1),1:n.x2] <- aperm(apply(l1, c(2,3), function(x) quantile(x, probs=quantiles)), c(2,1,3))
    l2[nrow(l2),1,1:n.x2] <- l2[n.x1,1,1:n.x2] - l2[1,1,1:n.x2]
    l2[nrow(l2),2:(n.q+1),1:n.x2] <- apply(l1[ ,n.x1,1:n.x2] - l1[ ,1,1:n.x2], 2, function(x) quantile(x, probs=quantiles))
    dimnames(l2) <- list(paste(c(rep(paste(x1name,"="),n.x1),paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),c(x1vals,"")),
                         c("est",quantiles),
                         paste(c(rep(paste(x2name,"="),n.x2)),
                               c(x2vals)))
    
    if (is.null(did)){did <- c(x2vals[1],x2vals[n.x2])} else{did <- did}
    l3 <- array(NA, c(1,n.q+1))
    l3[1,1] <- (l2[n.x1,1,match(did[2],x2vals)] - l2[1,1,match(did[2],x2vals)]) - (l2[n.x1,1,match(did[1],x2vals)] - l2[1,1,match(did[1],x2vals)])
    l3[1,2:(n.q+1)] <- quantile((l1[ ,n.x1,match(did[2],x2vals)] - l1[ ,1,match(did[2],x2vals)]) -  (l1[ ,n.x1,match(did[1],x2vals)] - l1[ ,1,match(did[1],x2vals)]), probs=quantiles)
    dimnames(l3) <- list("did",c("est",quantiles)) 
    
    ans <- new("post", 
               est=round(l2, digits=digits), 
               did=round(l3, digits=digits), 
               sims=l1, 
               model=class(model), 
               link=family(model)[2], 
               quantiles=quantiles, 
               call=call)
    return(ans)
  }
}


post.polr <- function(model,x1name=NULL,x1vals=NULL,x2name=NULL,x2vals=NULL,holds=NULL,seed=NULL,
                      n.sims=1000,cut=NULL,quantiles=c(.025,.975),did=NULL,weights=NULL, digits=3){
  
  call <- match.call()
  
  set.seed(seed)
  sims <- suppressMessages(postSim(model, n.sims=n.sims))
  
  if (is.null(weights)){wi <- c(rep(1, length(model$model[,1])))} else{wi <- weights}
  
  if (model$method=="probit"){link <- pnorm}
  else if (model$method=="logistic"){link <- plogis}
  else if (model$method=="cloglog"){link <- function(x){1-exp(-exp(x))}}
  else {stop("Link function is not supported")}
  
  n.obs <- length(model$model[,1])
  k <- length(model.matrix(polr(getCall(model)$formula, model$model))[1,])
  n.q <- length(quantiles)
  n.y <- length(levels(model$model[,1]))
  n.z <- length(model$zeta)
  tau <- array(NA, c(n.sims,n.z+2))
  zeta <- array(NA, c(n.z+2))
  tau[,1] <- zeta[1] <- -Inf
  tau[,2:(ncol(tau)-1)] <- sims@zeta[,1:n.z]
  zeta[2:(ncol(tau)-1)] <- model$zeta[1:n.z]
  tau[,ncol(tau)] <- zeta[ncol(tau)] <- Inf
  beta <- sims@coef
  
  if (is.null(cut)){
    
    if (is.null(x1name)){
      
      X_temp <- array(NA, c(n.obs,k))
      X <- array(NA, c(n.obs,k-1))
      
      newdata <- data.frame(model$model)
      if (!is.null(holds)){
        for (j in 1:length(holds)){
          newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
        }
      }
      X_temp[ , ] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
      X[ , ] <- X_temp[,-1]
      X <- aperm(X)
      
      l1 <- array(NA, c(n.sims, n.obs, n.y))
      l2 <- array(NA, c(n.sims,n.y))
      l3 <- array(NA, c(n.y,n.q+1))
      for (z in 1:n.y){
        l1[,,z] <- link(tau[,z+1] - beta %*% X) - link(tau[,z] - beta %*% X)
        l2[,z] <- apply(l1[,,z], 1, function(x) weighted.mean(x, wi))
        l3[z,1] <- weighted.mean(link(zeta[z+1] - coef(model) %*% X) - link(zeta[z] - coef(model) %*% X), wi)
        l3[z,2:(n.q+1)] <- quantile(l2[,z], probs=quantiles)
      }
      
      rownames(l3) <- paste(c(rep("Y =",n.y)), c(1:n.y))
      colnames(l3) <- c("est",quantiles)
      
      ans <- new("post", 
                 est=round(l3, digits=digits), 
                 did=NULL, 
                 sims=l2, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
    
    else if (is.null(x2name)){
      
      n.x1 <- length(x1vals)
      
      X_temp <- array(NA, c(n.obs,k,n.x1))
      X <- array(NA, c(n.obs,k-1,n.x1))
      
      for (i in 1:(n.x1)){
        newdata <- data.frame(model$model)
        if (!is.null(holds)){
          for (j in 1:length(holds)){
            newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
          }
        }
        newdata[ ,x1name] <- x1vals[i]
        X_temp[ , ,i] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
        X[ , ,i] <- X_temp[,-1,i]
      }
      
      l1 <- array(NA, c(n.sims, n.obs, n.x1, n.y))
      l2 <- array(NA, c(n.sims, n.x1, n.y))
      l3 <- array(NA, c(n.x1+1, n.q+1, n.y))
      X <- aperm(X, c(2,1,3))
      for (z in 1:n.y){
        l1[,,,z] <- apply(X, c(2,3), function(x) (link(tau[,z+1] - beta %*% x) - link(tau[,z] - beta %*% x)))
        l2[,,z] <- apply(l1[,,,z], c(1,3), function(x) weighted.mean(x, wi))
        l3[1:n.x1,1,z] <- apply(apply(X, c(2,3), function(x) (link(zeta[z+1] - coef(model) %*% x) - link(zeta[z] - coef(model) %*% x))), 2, function(x) weighted.mean(x, wi))
        l3[1:n.x1,2:(n.q+1),z] <- t(apply(l2[,,z], 2, function(x) quantile(x, probs=quantiles)))
        l3[nrow(l3),1,z] <- l3[n.x1,1,z] - l3[1,1,z]
        l3[nrow(l3),2:(n.q+1),z] <- quantile(l2[ ,n.x1,z] - l2[ ,1,z], probs=quantiles)
      }
      
      dimnames(l3) <- list(paste(c(rep(paste(x1name,"="),n.x1),paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),c(x1vals,"")),
                           c("est",quantiles),
                           paste(c(rep("Y =",length(levels(model$model[,1])))),
                                 c(1:length(levels(model$model[,1]))))) 
      
      ans <- new("post", 
                 est=round(l3, digits=digits), 
                 did=NULL, 
                 sims=l2, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
    
    else{
      
      n.x1 <- length(x1vals)
      n.x2 <- length(x2vals)
      
      X_temp <- array(NA, c(n.obs,k,n.x1,n.x2))
      X <- array(NA, c(n.obs,k-1,n.x1,n.x2))
      
      for (j in 1:n.x2){
        for (i in 1:(n.x1)){
          newdata <- data.frame(model$model)
          if (!is.null(holds)){
            for (k in 1:length(holds)){
              newdata[ ,names(holds)[k]] <- as.numeric(holds[k])
            }
          }
          newdata[ ,x1name] <- x1vals[i]
          newdata[ ,x2name] <- x2vals[j]
          X_temp[ , ,i,j] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
          X[ , ,i,j] <- X_temp[,-1,i,j]
        }
      }
      
      X <- aperm(X, c(2,1,3,4))
      
      l1 <- array(NA, c(n.sims, n.obs, n.x1, n.x2, n.y))
      l2 <- array(NA, c(n.sims, n.x1, n.x2, n.y))
      l3 <- array(NA, c(n.x1+1, n.q+1, n.x2, n.y))
      for (z in 1:n.y){
        l1[,,,,z] <- apply(X, c(2,3,4), function(x) (link(tau[,z+1] - beta %*% x) - link(tau[,z] - beta %*% x)))
        l2[,,,z] <- apply(l1[,,,,z], c(1,3,4), function(x) weighted.mean(x, wi))
        l3[1:n.x1,1,1:n.x2,z] <- apply(apply(X, c(2,3,4), function(x) (link(zeta[z+1] - coef(model) %*% x) - link(zeta[z] - coef(model) %*% x))), c(2,3), function(x) weighted.mean(x, wi))
        l3[1:n.x1,2:(n.q+1),1:n.x2,z] <- aperm(apply(l2[,,,z], c(2,3), function(x) quantile(x, probs=quantiles)), c(2,1,3))
        l3[nrow(l3),1,1:n.x2,z] <- l3[n.x1,1,1:n.x2,z] - l3[1,1,1:n.x2,z]
        l3[nrow(l3),2:(n.q+1),1:n.x2,z] <- apply(l2[,n.x1,,z] - l2[,1,,z], 2, function(x) quantile(x, probs=quantiles))
      }
      
      dimnames(l3) <- list(paste(c(rep(paste(x1name," ="),n.x1),paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),c(x1vals,"")),
                           c("est",quantiles),
                           paste(c(rep(paste(x2name,"="),n.x2)),x2vals),
                           paste(c(rep("Y =",n.y)), c(1:n.y))) 
      
      if (is.null(did)){did <- c(x2vals[1],x2vals[n.x2])} else{did <- did}
      l4 <- array(NA, c(n.y,n.q+1))
      for (i in 1:n.y){
        l4[i,1] <- (l3[n.x1,1,match(did[2],x2vals),i] - l3[1,1,match(did[2],x2vals),i]) -  (l3[n.x1,1,match(did[1],x2vals),i] - l3[1,1,match(did[1],x2vals),i])
        l4[i,2:(n.q+1)] <- quantile((l2[ ,n.x1,match(did[2],x2vals),i] - l2[ ,1,match(did[2],x2vals),i]) -  (l2[ ,n.x1,match(did[1],x2vals),i] - l2[ ,1,match(did[1],x2vals),i]), probs=quantiles)
      }
      yvals <- 1:n.y
      dimnames(l4) <- list(paste(c(rep(paste("Y","="),n.y)),yvals),c("est",quantiles)) 
      
      ans <- new("post", 
                 est=round(l3, digits=digits), 
                 did=round(l4, digits=digits), 
                 sims=l2, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
  }
  
  else{
    
    if (is.null(x1name)){
      
      X_temp <- array(NA, c(n.obs,k))
      X <- array(NA, c(n.obs,k-1))
      
      newdata <- data.frame(model$model)
      if (!is.null(holds)){
        for (j in 1:length(holds)){
          newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
        }
      }
      X_temp[ , ] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
      X[ , ] <- X_temp[,-1]
      X <- aperm(X)
      
      l1 <- apply(link(-tau[,cut+1] + beta %*% X), 1, function(x) weighted.mean(x, wi))
      l2 <- array(NA, c(1,n.q+1))
      l2[1,1] <- weighted.mean(link(-zeta[cut+1] + coef(model) %*% X))
      l2[1,2:(n.q+1)] <- quantile(l1, probs=quantiles)
      colnames(l2) <- c("est",quantiles)
      
      ans <- new("post", 
                 est=round(l2, digits=digits), 
                 did=NULL, 
                 sims=l1, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
    
    
    else if (is.null(x2name)){
      
      n.x1 <- length(x1vals)
      
      X_temp <- array(NA, c(n.obs,k,n.x1))
      X <- array(NA, c(n.obs,k-1,n.x1))
      
      for (i in 1:(n.x1)){
        newdata <- data.frame(model$model)
        if (!is.null(holds)){
          for (j in 1:length(holds)){
            newdata[ ,names(holds)[j]] <- as.numeric(holds[j])
          }
        }
        newdata[ ,x1name] <- x1vals[i]
        X_temp[ , ,i] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
        X[ , ,i] <- X_temp[,-1,i]
      }
      
      X <- aperm(X, c(2,1,3))
      l1 <- apply(apply(X, c(2,3), function(x) link(-tau[,cut+1] + beta %*% x)), c(1,3), function(x) weighted.mean(x, wi))
      l2 <- array(NA, c(n.x1+1,n.q+1))
      l2[1:n.x1,1] <- apply(apply(X, c(2,3), function(x) link(-zeta[cut+1] + coef(model) %*% x)), 2, function(x) weighted.mean(x,wi))
      l2[1:n.x1,2:(n.q+1)] <- t(apply(l1, 2, function(x) quantile(x, probs=quantiles)))
      l2[nrow(l2),1] <- l2[n.x1,1] - l2[1,1]
      l2[nrow(l2),2:(n.q+1)] <- quantile(l1[ ,ncol(l1)] - l1[ ,1], probs=quantiles)
      rownames(l2) <- c(paste(c(rep(paste(x1name,"="),n.x1),
                                paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),
                              c(x1vals,"")))  
      colnames(l2) <- c("est",quantiles)
      
      ans <- new("post", 
                 est=round(l2, digits=digits), 
                 did=NULL, 
                 sims=l1, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    } 
    
    else{
      
      n.x1 <- length(x1vals)
      n.x2 <- length(x2vals)
      
      X_temp <- array(NA, c(n.obs,k,n.x1,n.x2))
      X <- array(NA, c(n.obs,k-1,n.x1,n.x2))
      
      for (j in 1:n.x2){
        for (i in 1:n.x1){
          newdata <- data.frame(model$model)
          if (!is.null(holds)){
            for (k in 1:length(holds)){
              newdata[ ,names(holds)[k]] <- as.numeric(holds[k])
            }
          }
          newdata[ ,x1name] <- x1vals[i]
          newdata[ ,x2name] <- x2vals[j]
          X_temp[ , ,i,j] <- suppressWarnings(model.matrix(polr(getCall(model)$formula, data=newdata)))
          X[ , ,i,j] <- X_temp[,-1,i,j]
        }
      }
      
      X <- aperm(X, c(2,1,3,4))
      l1 <- apply(apply(X, c(2,3,4), function(x) link(-tau[,cut+1] + beta %*% x)), c(1,3,4), function(x) weighted.mean(x, wi))
      l2 <- array(NA, c(n.x1+1,n.q+1,n.x2))
      
      l2[1:n.x1,1,1:n.x2] <- apply(apply(X, c(2,3,4), function(x) link(-zeta[cut+1] + coef(model) %*% x)), c(2,3), function(x) weighted.mean(x, wi))
      l2[1:n.x1,2:(n.q+1),1:n.x2] <- aperm(apply(l1, c(2,3), function(x) quantile(x, probs=quantiles)), c(2,1,3))
      l2[nrow(l2),1,1:n.x2] <- l2[n.x1,1,1:n.x2] - l2[1,1,1:n.x2]
      l2[nrow(l2),2:(n.q+1),1:n.x2] <- apply(l1[ ,n.x1,1:n.x2] - l1[ ,1,1:n.x2], 2, function(x) quantile(x, probs=quantiles))
      
      dimnames(l2) <- list(paste(c(rep(paste(x1name," ="),n.x1),paste("\u0394","(",x1vals[1],",",x1vals[length(x1vals)],")")),c(x1vals,"")),
                           c("est",quantiles),
                           paste(c(rep(paste(x2name," ="),n.x2)),
                                 c(x2vals)))
      
      if (is.null(did)){did <- c(x2vals[1],x2vals[n.x2])} else{did <- did}
      l3 <- array(NA, c(1,n.q+1))
      l3[1,1] <- (l2[n.x1,1,match(did[2],x2vals)] - l2[1,1,match(did[2],x2vals)]) - (l2[n.x1,1,match(did[1],x2vals)] - l2[1,1,match(did[1],x2vals)])
      l3[1,2:(n.q+1)] <- quantile((l1[ ,n.x1,match(did[2],x2vals)] - l1[ ,1,match(did[2],x2vals)]) -  (l1[ ,n.x1,match(did[1],x2vals)] - l1[ ,1,match(did[1],x2vals)]), probs=quantiles)
      dimnames(l3) <- list("did",c("est",quantiles)) 
      
      ans <- new("post", 
                 est=round(l2, digits=digits), 
                 did=round(l3, digits=digits), 
                 sims=l1, 
                 model=class(model), 
                 link=model$method, 
                 quantiles=quantiles, 
                 call=call)
      return(ans)
    }
  }
}

setMethod("post", signature(model = "lm"), post.glm)
setMethod("post", signature(model = "glm"), post.glm)
setMethod("post", signature(model = "svyglm"), post.glm)  
setMethod("post", signature(model = "polr"), post.polr)  

    

